{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 parts of this notebook are from ([this Jupyter notebook](https://nbviewer.jupyter.org/github/krischer/seismo_live/blob/master/notebooks/Computational%20Seismology/The Finite-Difference Method/fd_first_derivative.ipynb)) by Kristina Garina and Heiner Igel ([@heinerigel](https://github.com/heinerigel)) which is a supplemenatry material to the book [Computational Seismology: A Practical Introduction](http://www.computational-seismology.org/),  additional modifications by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finite-difference approximations\n",
    "\n",
    "In the last lecture we covered different approaches to discretize continous media.\n",
    "Now, we only have to understand how to approximate (partial) derivatives by finite-differences."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Approximation of first derivatives\n",
    "\n",
    "The approximation of first derivatives using finite-differences is quite straightforward. Imagine that we have some time-dependent function $f(t)$, for which we want to calculate the first derivatives $\\frac{\\partial f(t)}{\\partial t}$:\n",
    "\n",
    "<img src=\"../images/gauss_disc_final.png\" width=\"95%\">\n",
    "\n",
    "First, we discretize $f(t)$ at discrete times\n",
    "\n",
    "\\begin{equation}\n",
    "t = i * dt, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "where $dt$ denotes a constant temporal sample interval. By definition, the first derivative of $f(t)$ with respect to $t$ is:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial f(t)}{\\partial t} = \\lim_{dt \\rightarrow 0} \\frac{f(t+dt)-f(t)}{dt}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "In the **finite-difference (FD) approximation**, we neglect the limit $dt \\rightarrow 0$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial f(t)}{\\partial t} \\approx \\frac{f(t+dt)-f(t)}{dt} \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Because $dt$ remains finite, this approximation is called finite-differences. Furthermore, we can distinguish between different finite difference operators, depending on the points involved in the approximation. The above example is a **forward FD operator**\n",
    "\n",
    "\\begin{equation}\n",
    "\\biggl(\\frac{\\partial f(t)}{\\partial t}\\biggr)^+ \\approx \\frac{f(t+dt)-f(t)}{dt}. \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "Alternatively, we can also define a **backward FD operator**\n",
    "\n",
    "\\begin{equation}\n",
    "\\biggl(\\frac{\\partial f(t)}{\\partial t}\\biggr)^- \\approx \\frac{f(t)-f(t-dt)}{dt}. \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "By taking the arithmetic average of the forward and backward operator\n",
    "\n",
    "\\begin{equation}\n",
    "\\biggl(\\frac{\\partial f(t)}{\\partial t}\\biggr) \\approx \\frac{\\left(\\frac{\\partial f(t)}{\\partial t}\\right)^- + \\left(\\frac{\\partial f(t)}{\\partial t}\\right)^+}{2},\\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "we get the **central FD operator**\n",
    "\n",
    "\\begin{equation}\n",
    "\\biggl(\\frac{\\partial f(t)}{\\partial t}\\biggr) \\approx \\frac{f(t+dt)-f(t-dt)}{2dt}. \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "##### Exercise\n",
    "\n",
    "Prove for the quadratic function\n",
    "\n",
    "\\begin{equation} \n",
    "g(t) = b t^2 \\nonumber\n",
    "\\end{equation} \n",
    "\n",
    "where $b$ is a constant time-independent parameter, that in the limit $dt \\rightarrow 0$, the forward, backward and central FD operators lead to the correct temporal first derivative:\n",
    "\n",
    "\\begin{equation} \n",
    "\\frac{\\partial g(t)}{\\partial t} = 2 b t \\nonumber\n",
    "\\end{equation} "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Import Libraries\n",
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, let's try to calculate the first temporal derivative of the Gaussian function\n",
    "\n",
    "\\begin{equation} \n",
    "f(t)=\\dfrac{1}{\\sqrt{2 \\pi a}}e^{-\\dfrac{(t-t_0)^2}{2a}},\\nonumber\n",
    "\\end{equation} \n",
    "\n",
    "where $a$ denotes the half-width of the Gaussian and $t_0$ a time-shift of the maximum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Define figure size\n",
    "rcParams['figure.figsize'] = 12, 5\n",
    "\n",
    "# Initial parameters\n",
    "   # maximum time\n",
    "   # number of time sample\n",
    "   # half-width      \n",
    "   # defining dt\n",
    "   # defining time shift t0\n",
    "\n",
    "   # defining time\n",
    "\n",
    "# Define gaussian function          \n",
    "\n",
    "    \n",
    "# Plotting of gaussian\n",
    "plt.plot(time, f)\n",
    "plt.title('Gaussian function')\n",
    "plt.xlabel('Time, s')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.xlim((0, tmax))\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Next, we calculate the numerical derivative using finite-differences with the forward operator:\n",
    "\n",
    "\\begin{equation}\n",
    "\\biggl(\\frac{\\partial f(t)}{\\partial t}\\biggr)^+ \\approx \\frac{f(t+dt)-f(t)}{dt}. \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "backward operator:\n",
    "\n",
    "\\begin{equation}\n",
    "\\biggl(\\frac{\\partial f(t)}{\\partial t}\\biggr)^- \\approx \\frac{f(t)-f(t-dt)}{dt}, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "and central operator:\n",
    "\n",
    "\\begin{equation}\n",
    "\\biggl(\\frac{\\partial f(t)}{\\partial t}\\biggr) \\approx \\frac{f(t+dt)-f(t-dt)}{2dt}. \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "To test the accuarcy of these approaches, we compare them with the analytical derivative:\n",
    "\n",
    "\\begin{equation} \n",
    "\\frac{\\partial f(t)}{\\partial t}=-\\dfrac{t-t_0}{a}\\dfrac{1}{\\sqrt{2\\pi a}}e^{-\\dfrac{(t-t_0)^2}{2a}} \\nonumber\n",
    "\\end{equation} "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# First derivative of Gaussian function\n",
    "\n",
    "# Initiation of numerical and analytical derivatives \n",
    "      # forward FD operator\n",
    "      # backward FD operator\n",
    "      # central FD operator\n",
    "      # analytical derivative\n",
    "\n",
    "# Numerical FD derivative of the Gaussian function\n",
    "for it in range (1, nt-1):\n",
    "    \n",
    "            # forward operator\n",
    "            # backward operator\n",
    "            # central operator\n",
    "\n",
    "# Analytical derivative of the Gaussian function     \n",
    "\n",
    "# Plot of the first derivative and analytical derivative\n",
    "plt.plot(time, nder_for,label=\"FD forward operator\", lw=2, color=\"y\")\n",
    "plt.plot(time, nder_back,label=\"FD backward operator\", lw=2, color=\"b\")\n",
    "plt.plot(time, nder_cent,label=\"FD central operator\", lw=2, color=\"g\")\n",
    "plt.plot(time, ader, label=\"Analytical derivative\", ls=\"--\",lw=2, color=\"red\")\n",
    "plt.title('First derivative')\n",
    "plt.xlabel('Time, s')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The approximation of the first derivative of the Gaussian by all three FD operators seem to be very accurate. How large are actually the errors?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Erros of the FD approximation\n",
    "\n",
    "# Plot of the errors of the first derivative FD approximations\n",
    "plt.plot(time, (nder_for-ader),label=\"FD forward operator error\", lw=2, color=\"y\")\n",
    "plt.plot(time, (nder_back-ader),label=\"FD backward operator error\", lw=2, color=\"b\")\n",
    "plt.plot(time, (nder_cent-ader),label=\"FD central operator error\", lw=2, color=\"g\")\n",
    "plt.plot(time, ader, label=\"Analytical derivative\",lw=2, color=\"red\")\n",
    "plt.title('First derivative errors')\n",
    "plt.xlabel('Time, s')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The errors for all three operators are actually quite small, so we have to amplify them by a factor 10. Larger errors occur in areas with significant slope. Notice also, that the central operator has the smallest error. You wonder why? We will see in a later lecture."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FD approximation of the 2nd derivative\n",
    "\n",
    "We derived FD approximations for the 1st derivative. However, to approximate the acoustic wave equation, we also need to know how to approximate a 2nd derivative. \n",
    "To achieve this, we calculate the **Taylor expansion** of $f(t+dt)$ up to the second order term:\n",
    "\n",
    "\\begin{equation}\n",
    "f(t+dt) \\approx f(t) + \\frac{\\partial f(t)}{\\partial t} dt + \\frac{1}{2}\\frac{\\partial^2 f(t)}{\\partial t^2} dt^2 + \\mathcal{O}(dt^3)\\nonumber\n",
    "\\end{equation} \n",
    "\n",
    "We do the same for $f(t-dt)$:\n",
    "\n",
    "\\begin{equation}\n",
    "f(t-dt) \\approx f(t) - \\frac{\\partial f(t)}{\\partial t} dt + \\frac{1}{2}\\frac{\\partial^2 f(t)}{\\partial t^2} dt^2 + \\mathcal{O}(dt^3)\\nonumber\n",
    "\\end{equation} \n",
    "\n",
    "Adding the approximations $f(t+dt)$ and $f(t-dt)$ leads to \n",
    "\n",
    "\\begin{equation}\n",
    "f(t-dt) + f(t+dt) \\approx 2 f(t) + \\frac{\\partial^2 f(t)}{\\partial t^2} dt^2\\nonumber\n",
    "\\end{equation} \n",
    "\n",
    "Finally, we can rearrange to $\\frac{\\partial^2 f(t)}{\\partial t^2}$ and get the following **FD approximation for the 2nd derivative**\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial^2 f(t)}{\\partial t^2} \\approx \\frac{f(t-dt) + f(t+dt) - 2 f(t)}{dt^2}\\nonumber\n",
    "\\end{equation} \n",
    "\n",
    "##### Excercise\n",
    "\n",
    "* Calculate the 2nd derivative $\\frac{\\partial^2 f(t)}{\\partial t^2}$ of the Gaussian \n",
    "\\begin{equation} \n",
    "f(t)=\\dfrac{1}{\\sqrt{2 \\pi a}}e^{-\\dfrac{(t-t_0)^2}{2a}},\\nonumber\n",
    "\\end{equation}\n",
    "analytically\n",
    "* Compute and compare the numerical and analytical 2nd derivative of the Gaussian and the errors between both solutions\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# 2nd derivative of Gaussian function\n",
    "\n",
    "# Initiation of numerical and analytical derivatives \n",
    "nder2=np.zeros(nt)      # 2nd derivative FD operator\n",
    "ader2=np.zeros(nt)      # analytical 2nd derivative\n",
    "\n",
    "# Numerical FD derivative of the Gaussian function\n",
    "for it in range (1, nt-1):\n",
    "    \n",
    "    # ADD YOUR 2ND DERIVATIVE OF THE GAUSSIAN HERE!\n",
    "    # nder2[it]=        # 2nd derivative FD operator\n",
    "    \n",
    "# ADD ANALYTICAL 2ND DERIVATIVE OF THE GAUSSIAN HERE!\n",
    "# ader2=\n",
    "\n",
    "# Plot of the numerical and analytical second derivative of the Gaussian\n",
    "#plt.plot(time, nder2,label=\"FD operator\", lw=2, color=\"y\")\n",
    "#plt.plot(time, ader2, label=\"Analytical 2nd derivative\", ls=\"--\",lw=2, color=\"red\")\n",
    "plt.title('2nd derivative Gaussian')\n",
    "plt.xlabel('Time, s')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "# Error of the 2nd derivative FD approximation\n",
    "\n",
    "# Plot of the second derivative FD approximation error\n",
    "#plt.plot(time, nder2-ader2,label=\"FD operator error\", lw=2, color=\"violet\")\n",
    "#plt.plot(time, ader2,label=\"Analytical 2nd derivative\", lw=2, color=\"red\")\n",
    "plt.title('2nd derivative errors')\n",
    "plt.xlabel('Time, s')\n",
    "plt.ylabel('Amplitude')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We learned:\n",
    "\n",
    "* Approximation of numerical first derivative of a time dependent function by forward, backward and central FD operators\n",
    "* Comparison with analytical solution to estimate the numerical errors \n",
    "* Approximation of 2nd derivative of a time dependent function using Taylor series expansion"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
